Quasi-equilibrium binary black hole sequences for 
puncture data derived from helical Killing vector conditions 
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We construct a sequence of binary black hole puncture data derived under the assumptions (i) 
that the ADM mass of each puncture as measured in the asymptotically flat space at the puncture 
stays constant along the sequence, and (ii) that the orbits along the sequence are quasi-circular in 
the sense that several necessary conditions for the existence of a helical Killing vector are satisfied. 
These conditions are equality of ADM and Komar mass at infinity and equality of the ADM and a 
rescaled Komar mass at each puncture. In this paper we explicitly give results for the case of an 
equal mass black hole binary without spin, but our approach can also be applied in the general case. 
We find that up to numerical accuracy the apparent horizon mass also remains constant along the 
sequence and that the prediction for the innermost stable circular orbit is similar to what has been 
found with the effective potential method. 

PACS numbers: 04.20.Ex 04.25. Dm, 04.30.Db, 95.30.Sf 
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I. INTRODUCTION 

Binary black hole inspirals and mergers are promising 
sources for ground-based interferometric gravitational 
wave detectors such as GEO600, LIGO and TAMA [ij. 
These systems are highly relativistic once they enter the 
sensitive frequency band of the detector. Hence, in order 
to predict the gravitational waves from these sources nu- 
merical simulations will be needed, which in turn require 
astrophysically realistic initial data. Several methods to 
produce initial data for binary black hole systems exist 
(see for a recent review). It is, however, not clear yet 
which method should be used to obtain initial data that 
are as realistic as possible in terms of astrophysical con- 
tent. During the inspiral, Post-Newtonian (PN) theory 
predicts that the two black holes will be in quasi-circular 
orbits around each other with a radius which shrinks on 
a timescale much larger than the orbital timescale. One 
approach to incorporate this PN information is to di- 
rectly use the PN metric when solving the constraint 
equations to obtain initial data (sj. Another approach 
is the so called effective potential method which iden- 
tifies circular orbits with minima in a suitably defined 
binding energy Yet another approach is to concen- 
trate on the fact that quasi-circular inspiral orbits evolve 
slowly, which means that the initial data should have an 
approximate helical Killing vector |^ . 

Not only are there different methods to define quasi- 
circularity, but there are different ways to construct bi- 
nary black hole initial data on arbitrary orbits. Due to its 
simplicity, the puncture construction of black hole data 
and the puncture evolution method ^ have been used 
in many black hole simulations, in fact to date all gravita- 
tional wave forms obtained numerically for binary black 
hole inspirals arc based on puncture initial data, see e.g. 
^, 10, 11, 12, 13, 14]. Although there are several al- 
ternatives (21 J concentrate on puncture data in this 
paper. 



The first quasi-circular puncture sequence was ob- 
tained by Baumgarte using the effective potential 
method and assuming constant apparent horizon mass 
along the sequence. Baker derived a puncture se- 
quence based on a variational principle that leads to 
an effective potential method with constancy of the 
Arnowitt-Deser-Misner (ADM) mass at the puncture. 
Furthermore, within a certain error a puncture sequence 
can also be obtained |Q| by using the orbital parame- 
ters computed for Misner type black hole excision data 
by Cook g. 

Recently, we have shown how the approximate heli- 
cal Killing vector idea can be applied to puncture data 
iSi]. We have found that puncture data with orbital pa- 
rameters determined using the effective potential method 
approximately fulfill several necessary conditions for the 
existence of a helical Killing vector, namely the equality 
of the ADM and Komar mass at infinity, and equality of 
ADM mass and a suitably scaled Komar mass at both 
punctures. 

The main result of the present paper is that this work 
can be extended to obtain sequences of black hole punc- 
ture data in quasi-circular orbits. We use the same neces- 
sary conditions derived for an approximate helical Killing 
vector plus the assumption that the ADM mass of each 
puncture as measured in the asymptotically flat space at 
the puncture stays constant. Our results for the inner- 
most stable circular orbit (ISCO) are close to what P| 
and jl5|| find with the effective potential method. We 
also find that up to numerical accuracy the apparent 
horizon mass stays constant along our sequence, which 
is usually assumed to hold in conjunction with the effec- 
tive potential method. Furthermore, we provide detailed 
information for the orbital parameters of the puncture 
sequence, which has not yet been available in the litera- 
ture since 0,^^ focus on the ISCO and do not provide 
this information. 
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II. QUASI-CIRCULAR SEQUENCES 



The orbits of two puncture black holes are completely 1 .000 



masses normalized to M=1 



described by the ADM mass Mj^ and spin S\ of each 
black hole computed in the asymptotically flat space at 
the puncture, the coordinate distance D between the two 
black holes, and the the momentum parameter of 
each black hole {A = 1,2 for the two black holes). If we 
choose coordinates such that the total ADM momentum 
is zero, the two black holes must have equal and oppo- 
site momenta, and hence only one momentum parameter 
P = \P\\ is needed. Thus the system is characterized by 
the 10 physical parameters M^^*^, S\, P, and D. In ad- 
dition, we have the 2 parameters ca which fix the value 
of the lapse a at each puncture. In our approach a is 
needed to compute Komar mass integrals at both punc- 
tures and at infinity, and we determine a from a maximal 
slicing equation Hence, puncture data together with 
the lapse are described by 12 parameters. 

We want to construct a sequence of quasi-circular or- 



bits, 



i.e. we want to find M^"^ , 



S\, P, and Ca such 



that for any D we have the same two black holes on a 
quasi-circular orbit. In order to assure that we are always 
dealing with the same two black holes, we impose that 
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FIG. 1: The ADM mass at infinity, the total 2PN energy, and 
the apparent horizon mass as a function of angular velocity. 
The minimum in M^^^' at Mn^in = 0.19 ± 0.03 is marked 
by a dot. 



S\ — const. 



(2) 



Note that this choice is not unique, since individual 
masses and spins of black holes in binaries contain am- 
biguities. But it is expected that Mj^^*^ and Sa will 
be approximately constant during the inspiral of realis- 
tic black holes, at least in the quasi-stationary regime 
[TtI. To assume constancy of the masses at the punc- 
tures was also considered by Baker 0. In addition, to 
ensure quasi-circular orbits, we impose the three neces- 
sary helical Killing vector conditions |^ 



ADM 



and 



Here 



Ik (a, Sa) 



lK{a,S) 



„ n.rADM 
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(3) 



(4) 
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is the Komar mass, evaluated either on a sphere Soo at 
infinity or on infinitesimal spheres Sa around each punc- 
ture. Now the 11 parameters M^^^^ , S^^, P, and ca can 
be determined from the 11 conditions Q, ©, and 
Q for any given black hole separation D. 

For simplicity, we now restrict ourselves to non- 
spinning equal mass binaries with M^^'^'^ — AI/2. In 
this case ci = C2 = c will hold, so that we are left with 
only the 3 parameters M , P, and c, which can be fixed 



using the 3 conditions Q, and Q. The number of 
parameters can be further reduced through rescaling by 
M. When this is done it is sufficient to find the 2 pa- 
rameters P/AI and c which satisfy conditions Q and Q 
for each given D/M. Determining the appropriate pa- 
rameters involves root finding, which is computationally 
expensive since in each iteration two elliptic equations 
for the conformal factor (j) and for the lapse a have to be 
solved. 

In order to construct sequences of quasi-circular orbits 
as described above, we use the same numerical techniques 
as in combined with a Newton-Raphson method to 
find the parameters P/M and c. All numerical grids 
have a uniform resolution of h = 0.0625M, and the outer 
boundary is located at Ri — 12M. To obtain accu- 
rate values for the ADM and Komar integrals needed 
in Eq. (PJ we use volume integrals over the numerical 
grid (out to Ri — 12M) plus a correction from inte- 
grating outside the grid out to a larger distance i?2 (see 
appendix ^) . To sufficiently reduce the error with this 
correction we set i?2 = lORj/M = 1440M. Then the 
masses at infinity have errors of only 0.0005%. The ADM 
and Komar integrals at the punctures needed in condi- 
tion Q are obtained by fourth order interpolation of the 
conformal factor and the lapse onto the location of the 
punctures. The error in the masses at the punctures is of 
order 0.03%, which is thus also the error to within which 
we can determine our sequence, as it is much larger than 
the error in the masses at infinity. 



In Fig.nwe show our results for the ADM mass M, 



ADM 
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at infinity and the apparent horizon mass 



(6) 



Here Ai and A2 are the areas of the two apparent hori- 
zons. In order to obtain coordinate independent results, 
all our graphs are plotted versus the angular velocity ^ 

O _ M^D^' ~lK{a,S,)~lK{o^,S2) 

"~ 2JADM ' 

instead of the coordinate separation D. For comparison, 
we also show the 2PN total energy 0, IT^ 



■2PN 



M J AD 
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versus the 2PN angular velocity 0, [l8j | 
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M(D/Mf 
(1 + 2D/Mf 
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Note that VI2PN in Eq- is written such that ri2PAf is 
exact up to all PN orders in the limit of /i/A/ ^ 0, while 
for 11 /M > Eq. © is accurate up to 2PN order. It 
should be kept in mind, however, that both, the resum- 
mation of Eq. @ , as well as the fact that E2P n is given 
as a function of D instead of say fl2PN, will change the 
E2PN versus VI curve by PN terms higher than the 2PN 
terms included. This effect becomes noticeable around 
MQ, > 0.1, where these higher order terms become im- 
portant. 

The ADM mass Af^^*^ approaches the PN result for 
large separations (i.e. small Vl) as expected. Also, within 
our accuracy M^^ is constant along the sequence. In 
previous works 0, IT^j l this property was usually enforced 
in place of Eq. It is interesting that the necessary 
conditions and Q for the existence of a helical Killing 
vector imply that Af^^*^ = const and M^^ = const are 
equivalent, and that in fact 2Af^^*^ = Af = M^" . The 
ADM mass M^°^ has a minimum at Mfl = Mn^i„ = 
0.19 ±0.03, which is usually interpreted as the innermost 
stable circular orbit (ISCO). This ISCO minimum occurs 
(within error bars) at the same angular velocity as the 
ones found by Cook 3 ^nd Baumgarte using the ef- 
fective potential method. Hence, results obtained using 
the assumption of a helical Killing vector are very simi- 
lar to results found using the effective potential method. 
The reason for the relatively large error bar in Mflmin 
comes from the fact that the minimum in M^'^ is very 
shallow, so that the errors of order 0.03% in our mass de- 
termination get amplified. Note, however, that our error 
estimate for MVLmin is very conservative, as we have as- 
sumed that Af^^*^ might randomly vary within 0.03%. 
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FIG. 2: The ADM angular momentum at infinity and the 
2PN angular momentum as a function of angular velocity. 
The minimum in J^^" at Mf2™^i„ = 0.190 ±0.015 is marked 
by a dot. 



If we instead assume that all values for AT^^*^ are sys- 
tematically too large or too small by 0.03%, the errorbar 
would be much smaller. 

In Fig. 12 we show the ADM angular momentum J^^'^^ 
as well as the 2PN angular momentum 0, E| 
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versus Alil, which agree for small fl as expected. As 
we can see, has its minimum at the same f2 as 

Af^"^^^. However, as J^^^^ varies more strongly than 
M^^^ the minimum is less shallow, which allows us to 
locate it with better accuracy, with the result AfJ7 = 
MVtmin = 0.190±0.015. Note that even though the error 
in Mflrnin is relatively large, we are able to determine all 
other quantities with an error of only about 0.03%. 

Next, we list fits of several quantities along our numer- 
ically computed sequence. All fitted quantities are given 
as functions of the dimensionless orbital radius 



D 

2M' 



(11) 



for the case of an equal mass binary. We find 

TO 1 0.0668549 0.0280090 0.0906937 

M 



2 X xi 

0.139304 0.0973014 0.0238708 



(12) 



a;2 X 
P _ 1 0.0123341 0.274250 0.317647 

M ~ 4^/2iP, X 



x^ 
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0.0243978 



0.509913 



IX 



1.00585 



X 
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For the given expansions the fitting error is about 5 x 
10~^. Our fits can be used to put two punctures into a 
quasi-circular orbit for any desired separation D. For the 
reader's convenience we also list the results of these fits 
for a variety of separations in Tab. ^ 



III. DISCUSSION 

We have constructed a sequence of binary black hole 
puncture data by using the necessary helical Killing vec- 
tor conditions Q and Q and by assuming that M^^^^ 
is constant along the sequence. The numerical results 
are obtained for equal mass binaries without spin, but 
our approach can also be applied in the general case. 
Our sequence is close to quasi-equilibrium in the follow- 
ing sense. The time derivative of the trace of the extrinsic 
curvature is zero, due to the choice of a maximal slicing 
lapse. The conformal metric and conformal factor evolve 
on a timescale long compared to the orbital timescale, 
if we also compute a shift as in 0, which removes the 
longitudinal piece of the time derivative of the conformal 
metric. With this gauge choice the time derivative of the 
tracefree part of the extrinsic curvature is also reduced, 
but it still evolves on the orbital timescale. 

We find that the apparent horizon mass Mah is con- 
stant along the sequence up to our numerical accuracy, 
and that Mah = M^°^ + M^°^ . Furthermore, the 
ISCO minimum in both M^^^ and J^^*^ occurs at 
the same 57, but the error bars on the minimum on in 
J^^^^ are tighter. Previous results for the ISCO H ^ 
agree with our result within error bars, but our ISCO is 
at a somewhat higher fi. 

Baker has constructed a puncture sequence using 
an effective potential method and the assumption that 



j\^ADM jg constant along the sequence. His results are 
similar to ours. Yet even though he has used fixed mesh 
refinement for computational efficiency, his values seem 



TABLE I: Parameters for quasi-circular orbits of binary black 
hole puncture data. The errors in the quantities listed are of 
order 0.03%. 



D/2M 


m/M 


P/M 


c 




jAUM 
''oo 


MQ. 


M 




1.0300 


0.43193 


0.35879 


0.69584 


0.97681 


2.9564 


0.19008 


1.0500 


0.43318 


0.35196 


0.70109 


0.97682 


2.9564 


0.18653 


1.0750 


0.43469 


0.34380 


0.70743 


0.97683 


2.9567 


0.18223 


1.1000 


0.43613 


0.33604 


0.71354 


0.97685 


2.9571 


0.17809 


1.1250 


0.43752 


0.32865 


0.71942 


0.97689 


2.9578 


0.17409 


1.1500 


0.43885 


0.32160 


0.72509 


0.97693 


2.9587 


0.17024 


1.1750 


0.44012 


0.31488 


0.73055 


0.97698 


2.9599 


0.16652 


1.2000 


0.44135 


0.30846 


0.73581 


0.97703 


2.9612 


0.16292 


1.2250 


0.44253 


0.30232 


0.74089 


0.97710 


2.9627 


0.15945 


1.2500 


0.44366 


0.29645 


0.74580 


0.97717 


2.9645 


0.15610 


1.3000 


0.44581 


0.28543 


0.75511 


0.97732 


2.9684 


0.14972 


1.3500 


0.44780 


0.27529 


0.76380 


0.97750 


2.9731 


0.14375 


1.4000 


0.44965 


0.26594 


0.77194 


0.97769 


2.9785 


0.13815 


1.4500 


0.45138 


0.25728 


0.77957 


0.97789 


2.9844 


0.13290 


1.5000 


0.45300 


0.24924 


0.78673 


0.97811 


2.9909 


0.12797 


1.6000 


0.45593 


0.23480 


0.79981 


0.97856 


3.0055 


0.11895 


1.7000 


0.45853 


0.22219 


0.81144 


0.97903 


3.0217 


0.11092 


1.8000 


0.46084 


0.21108 


0.82184 


0.97951 


3.0395 


0.10375 


1.9000 


0.46291 


0.20123 


0.83120 


0.98000 


3.0586 


0.097299 


2.0000 


0.46477 


0.19243 


0.83964 


0.98048 


3.0788 


0.091487 


2.2500 


0.46870 


0.17406 


0.85755 


0.98164 


3.1330 


0.079215 


2.5000 


0.47185 


0.15956 


0.87192 


0.98272 


3.1911 


0.069443 


2.7500 


0.47442 


0.14780 


0.88369 


0.98372 


3.2517 


0.061518 


3.0000 


0.47656 


0.13808 


0.89350 


0.98461 


3.3138 


0.054988 


3.5000 


0.47992 


0.12287 


0.90891 


0.98619 


3.4404 


0.044923 


4.0000 


0.48243 


0.11148 


0.92044 


0.98750 


3.5674 


0.037589 


4.5000 


0.48439 


0.10259 


0.92940 


0.98859 


3.6934 


0.032052 


5.0000 


0.48595 


0.095433 


0.93655 


0.98952 


3.8173 


0.027753 


6.0000 


0.48830 


0.084541 


0.94726 


0.99099 


4.0580 


0.021565 


7.0000 


0.48997 


0.076578 


0.95491 


0.99211 


4.2884 


0.017378 


8.0000 


0.49123 


0.070451 


0.96064 


0.99299 


4.5089 


0.014390 


9.0000 


0.49220 


0.065559 


0.96511 


0.99369 


4.7203 


0.012171 


10.000 


0.49299 


0.061542 


0.96868 


0.99427 


4.9233 


0.010468 



to be less accurate. In the outer boundary was put 
at large distances in order to compute M^^'^ with suf- 
ficient accuracy, but based on our numerical results the 
resolution near the black holes appears to be rather low. 

The parameters for quasi-circular orbits in Tab. Q] can 
be compared with the data for the effective potential 
method in the form given in |l4j . which provides a se- 
quence translated from Misner type black hole excision 
data of Cook 0, to puncture data assuming that 
both types of data are numerically close. Our data is 
likely to be more accurate than effective potential method 
data since we do not have to find turning points in po- 
tentially very flat ADM mass curves to determine if an 
orbit is circular. Also, our masses are computed with 
much higher accuracy than in For example, for 
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D/2 = 1.7846M = 1.8490M^^*^ there is a discrepancy 
of 1%. Yet part of the difference may be due to the fact 
that we use the hehcal Kilhng vector conditions © and 
Q instead of effective potential method to find quasi- 
circular orbits. The fact that we have kept Aff con- 



stant along the sequence and not Mah as in |J| and 
should not play a role as both seem to be equivalent along 
our sequence. 

Finally, we want to mention that G. Cook has informed 
us that M. Hannam has just completed his Ph.D. thesis 
which includes the discussion of a puncture sequence that 
shares some of the features of our work. 
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APPENDIX A: COMPUTING MASSES AT 
INFINITY ON A FINITE GRID 

In order to implement Eqs. Q and I^J, we have to 
determine the masses with sufficient accuracy. In partic- 
ular, we have to compute the ADM mass at infinity to a 
rather high accuracy of about lO-'^M. Since in the case 
of punctures the 3-metric is conformally fiat with a 
conformal factor 



= 1 



mi m2 

the ADM mass at infinity is given by 
1 

2^ 



(Al) 



M, 



ADM 



1 

2^ 



m 

Here u has to be determined from an elliptic equation of 
the form 0,0 



with boundary condition 



lim u ^ 0. 



(A3) 



(A4) 



As our numerical grid does not extend all the way to in- 
finity, we cannot directly compute the volume integral in 
Eq. (|A2ll . Using a coordinate transformation in the ra- 
dial direction, e.g. of the "FishEye" type helps but 
does not necessarily give sufficiently accurate results be- 
cause of finite difference errors in the far region. Instead, 
the last term of Eq. (jA2|l can be approximated by the 
integral 



fiu) dV 



(A5) 
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FIG. 3: M^^*^ computed using integral ljA6|l over the nu- 
merical grid only, compared with M^^'^' computed using the 
extended volume integral 1 A 1411 with R2 = \QR\/M for two 
different resolutions h. The extended integral 1 A 1411 gives 
much more accurate results. 

over the numerical grid up to the outer boundary denoted 
by Ri. Then 



V^u dV = /i +0(l/i?i). 



(A6) 



where 0(l/i?i) is the error due to truncating the integral 
at a finite radius Ri. This error is often too large as 
it falls off only like \/Ri (see the curves labeled "grid 
only" in Fig.|2Jl. For this reason we want to improve the 
approximation by integrating over a region larger than 
the grid, which extends to at least R2 ^ Rf, so that the 
error goes down by one additional power in l/Ri- I.e. we 
want to approximate the volume integral by 



V^u dV^Ii+l2 + 0(1/ R2) ^h+h 



where 



I2 = I fiu) dV. 



0(l/i??), 
(A7) 

(A8) 



The problem with this approach is that we have com- 
puted u only on the numerical grid, so that we have to 
somehow approximate u outside the grid. This can be 
achieved by noting that because of boundary condition 
[1 



Ol 

r 



0{l/r^) 



(A9) 
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for large distances r, so that 

Ii = [ V^u dV ^ (f V,M dS' = -Anai + 0(l/i?i). 

Jq JR.I 

(AlO) 

Hence outside the numerical grid u can be approximated 
by 

where Ji is computed only on the grid. Also, note that 
for large r 

!{u) ^ ^, (A12) 

otherwise the ADM mass integral in Eq. (|A2|) would be 
infinite. Thus, if we insert Eq. I|A11(I into Eq. (|A8|) and 
also take into account the fallofF behavior ljA12p of f{u) 
we arrive at 



Combining Eqs. and (jA13|l we arrive at the final 

result 

/v=„.,.^7,+/*/(-±i).^ + 0(i,), 

(A14) 

which depends only on Ii computed from u on the nu- 
merical grid and a volume integral outside the grid over 
the known function /. 

We have applied the method described in the previ- 
ous paragraph to both the ADM and Komar mass inte- 
grals at infinity, as both can be expressed as volume inte- 
grals. Fig.|2|shows M^-^*^ computed using the extended 
volume integral (|A14|I as well as the result if we inte- 
grate only over the grid using Eq. ljA6p . The extended 
volume integral result is clearly more accurate, say for 
Ri = 12M, which we have used to obtain our sequence. 
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